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Here we study a driven lattice gas model for microtubule depolymerizing molecular motors, where 
traffic jams of motors induce stochastic switching between microtubule growth and shrinkage. We 
term this phenomenon “traffic dynamic instabiiity” because it is reminiscent of microtubule dynamic 
instability [T. Mitchison and M. Kirschner, Nature 312, 237 (1984)]. The intermittent dynamics of 
growth and shrinking emerges from the interplay between the arrival of motors at the microtubule 
tip, motor induced depolymerization, and motor detachment from the tip. The switching dynamics 
correlates with low and high motor density on the lattice. This leads to an effectively bistable 
particle density in the system. A refined domain waif theory predicts this transient appearance of 
different phases in the system. The theoretical results are supported by stochastic simulations. 


I. INTRODUCTION 


Microtubule (MT) depolymerizing enzymes [IH3] are 
considered important for MT length-regulation @HS]. 
These enzymes function in parallel to MT dynamic insta¬ 
bility [7], which is the hydrolysis-driven stochastic switch¬ 
ing of MTs between a growing and a shrinking state Eng. 
The class of MT depolymerizing enzymes [3] contains the 
kinesin-8 protein family, which are molecular motors that 
walk towards the MT plus-end [T0l - fl2l . These have been 
studied in detail in the biological literature, see Ref. [T3 
and references therein. 

Suitable theoretical models to describe the collective 
movement of molecular motors on filaments are driven 
lattice gases □MI]- Such models explain the formation 
of traffic jams on MTs [HMH] as observed in experi¬ 
ments [2T, 22j. Recently also lattice gases of dynamic 
system size were studied. For example when individual 
particles trigger lattice growth by forming new lattice 
sites at the lattice end [23IEfi] , or remove lattice sites HU¬ 
ES]. 

The interplay between lattice growth and shrinking 
was investigated in a variety of settings, e.g. when par¬ 
ticles stabilize shrinking lattices ES9. or depolymerize 
growing lattices [-T* [61, : 3JJ, [32]. An analogy to these molec¬ 
ular motor systems can be found in queuing theory, where 
the length of a queue is dynamic [33] and waiting times 
are of central interest [33]. 

Here we study a simplified stochastic model for a dy¬ 
namic MT tip and MT depolymerizing molecular mo¬ 
tors. The model is based on the totally asymmetric 
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simple exclusion process (TASEP) [35j and includes spe¬ 
cific stochastic processes that account for the growth and 
shrinking of the lattice. Particularly, we study a kinetic 
model where the following MT tip related processes are 
included: 1) motors depolymerize the MT at the tip, 2) 
motors inhibit MT growth if bound to the tip, and 3) 
motors have a finite dwell time at the tip. Remarkably, 
the system investigated here displays stochastic switch¬ 
ing between extended periods of MT growth and MT 
shrinking. This intermittent dynamics is reminiscent of 
MT dynamic instability [5] , yet its origin is different and 
relies on the dynamics of motors at the MT tip: We iden¬ 
tify the spontaneous formation of motor traffic jams at 
the MT tip as a “catastrophe” event, because it leads to 
shrinking of the MT. On the other hand, we identify the 
stochastic detachment of a motor from the tip as a “res¬ 
cue” event, because it initiates MT growth in the model. 
The predictions of the model can be tested in biochemical 
reconstitution experiments [36] and could help to identify 
detailed MT-motor interactions as they become increas¬ 
ingly accessible to experiments, see e.g. Ref. [ 37] , 


This work is organized as follows. In section [TT] the 
details of the model are presented. In the results section 
we present the phenomenon of “traffic dynamic instabil¬ 
ity” (III A). Then, the mean field solution of the model 
is introduced and discussed (III B). In sections III C and 


HID we numerically study the formation of shocks and 


develop a domain wall theory which quantifies our obser¬ 
vations. In particular, this domain wall theory allows to 
identify metastable regimes in the phase space of the sys¬ 


tem as is shown in the following sections: In section III E 


we identify a phase of stripe formation due to particu¬ 
lar velocities of domain wall motion in the system. The 
intermittent phases are analyzed in greater detail analyt¬ 


ically as well as numerically in sections IIIF and |III G 
Finally we discuss our results and conclude in section |IV 
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FIG. 1. Molecular motors on the MT lattice. At the left 
motors enter the lattice from a constant reservoir p_, which 
mimics the motor density in bulk of the MT. The motor den¬ 
sity at the tip (MT plus-end) is denoted p+. Motors detach 
from the MT tip at rate /3. The MT grows in absence of a 
motor at rate p and is depolymerized by a motor at rate 5. 
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II. MODEL 

We consider a lattice gas model for molecular motors 
close to the MT plus-end as illustrated in Fig. [T] Mo¬ 
tors move from the bulk of the MT (left) to the plus-end 
(right) at rate v = 1 if the next site is empty. This choice 
of v sets the timescale of all other rate constants. The 
dynamics can be formulated in terms of occupation num¬ 
bers Hi , where n % = 1 and rq = 0 denote the presence and 
absence of a particle, respectively. At the left, the sys¬ 
tem is coupled to a constant particle reservoir p_, which 
emulates the bulk of the MT. At the right, the model 
corresponds to the MT plus-end, where particles detach 
at rate /3. The above model is known as the totally asym¬ 
metric simple exclusion process (TASEP) and plays a 
paradigmatic role in non-equilibrium statistical mechan¬ 
ics \mm\ . Furthermore, we denote the probability that 
a motor occupies the tip-site with p+ and the probability 
to find a motor in bulk of the lattice pb ■ We consider a 
lattice of constant size N, which is co-moving with the 
plus-end of the MT. At the left boundary, the system 
looses a particle if the lattice grows and the leftmost site 
is occupied. In case of a depolymerization event, the lat¬ 
tice moves to the left, where one lattice site is added. The 
newly added site is occupied with probability p _. At the 
right boundary, two more processes can happen: (i) Par¬ 
ticles remove the terminal site with rate S [22J. fa) If the 
tip-site is empty the lattice grows at rate p [6L These 
dependencies of MT growth and shrinking on motor oc¬ 
cupation at the tip can also be interpreted as the two 
different nucleotide states of MT tips [5]: in the GTP 
state the filament grows, corresponding to the absence 
of a motor in our model; in the GDP state the filament 
shrinks, corresponding to a motor bound to the tip in 
our model. The dynamics of motors on the lattice, how¬ 
ever, renders our model different from GTP hydrolysis 
dependent switching, and, therefore, distinct from MT 
dynamic instability. In the model presented here switch¬ 
ing between growth and shrinking (catastrophe) is due 
to the spontaneous formation of traffic jams on growing 
MTs. Switching from shrinking to growth (rescue) is due 
to detachment of a motor from the tip. 


FIG. 2. Stochastic simulations show intermittent dynamics. 
Switching between periods of growth and shrinking are rem¬ 
iniscent of MT dynamic instability and indicate bistable be¬ 
havior of the system: The three panels show the trajectory of 
the lattice tip, the bulk motor density on the MT pb, and the 
motor occupation of the tip n+. Shaded areas indicate peri¬ 
ods of growth that involve low bulk density pb ss 0.2. System 
size is N — 200. 


III. RESULTS 
A. Intermittent dynamics 

In stochastic simulations of the model we observe in¬ 
termittent dynamics. In Fig. [2] this is illustrated for 
three key observables, which are the trajectory of the 
MT tip, the average bulk density of motors on the lat¬ 
tice pb = n ii an( I Ih e motor occupation at the 

tip n+ € {0; 1}, see Fig. [2j Periods of depolymeriza¬ 
tion are distinguished from periods of growth by a high 
bulk density pb « 1 and persistent occupation of the 
tip ri+ = 1, cf. the white and shaded areas in Fig. 2] 
respectively. The underlying switching processes of the 
system can be understood intuitively in terms of a sepa¬ 
ration of timescales: If the tip site is occupied, n+ = 1, it 
constitutes a bottleneck [421 155] behind which particles 
pile up due to motor traffic and slow depolymerization. 
For rare tip detachment (/? -C 1) this bottleneck persists 
and eventually induces a high density of motors in bulk, 
Pb ~ 1. The high bulk density manifests the shrinking 
state of the filament. However, the system may stochasti¬ 
cally switch to growth. This can happen on the occasion 
that the bottleneck at the tip site is removed through tip 
detachment (rate /?). Then the tip site becomes empty, 
n+ = 0, which entails that the lattice is in its growing 
state. The timescale of this rescue event is /3 _1 . 
Similarly, it is possible to obtain a qualitative under¬ 
standing of catastrophe events: At high growth rates p, 
motors can not keep up with the growing tip. The system 
assumes a steady state where the bulk density is signifi¬ 
cantly smaller than in the depolymerizing state, pb ~ 0.2, 
see Fig. [2] Stochastic switching into the shrinking state 
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FIG. 3. (a) Comparison of stochastic simulations (symbols) 
and analytic results (lines) obtained for the tip densities from 
mean-held calculations; see Eqs. Q, (J5|, and |7|. For rj = 
0.35, p_ = 0.5, /3 = 0.05, and N = 400 we find very good 
agreement between theory and simulations. The transition 
between the IN and the MC phase can hardly be recognized; 
we plotted both functions and elaborate on the exact criterion 
later, (b) Data for bulk densities (symbols) also confirms the 
validity of mean-held calculations [lines as in panel (a)]. The 
transition between the EX and the IN phase is discontinuous 
in the bulk densities in agreement with what is known from 
TASEP: At the transition between the EX and the IN phase, 
both phases coexist and are separated by a diffusing domain- 
wall (DW). (c) The DW ensues density fluctuations, which we 
measured in terms of the normalized standard deviation of the 
bulk density. At the EX/IN-transition density huctuations 
show a pronounced peak (indicated by the vertical line). 


may happen if a single motor reaches the tip site and 
induces traffic jam formation, which reverses growth to 
shrinking. Thus, the timescale of catastrophe is given by 
the arrival rate of motors at the tip, j + . This arrival rate 
is sensitive to current fluctuations and the system size 

n mi. 


B. Limitations of the mean-field approach 


In a first step to understand the dynamics of the sys¬ 
tem, we determine the tip-densities, p+, within a mean- 
field approach [5j. There are three generic phases of 
motor dynamics on the lattice: a low density phase 
(IN phase), a high density phase (EX phase), and a 
maximal current phase (MC phase). The dynamics of 
such a system is likewise dependent on particle input 
p_, the particle exit rate /J, or the capacity of particle 
flow on the lattice, respectively l38lTTT| . As shown re¬ 
cently, this requires an analysis of bulk and boundary 
currents in the system [5]- Employing a mean-field ap¬ 
proximation for nearest neighbor occupation numbers, 
(riini+i) = (ni)(ni+i), the bulk current of the system 
reads 0 

MPb, P+) = Pb(l - Pb) + SpbP+ ~ ?7Pb(l - P+), (1) 

where the terms on the right hand side stand for par¬ 
ticle hopping with on-site exclusion, depolymerization, 
and polymerization of the MT. Note that the latter de¬ 
pend on the probability that a motor is bound to the 
tip. In terms of particle movements, depolymerization 
and polymerization correspond to parallel updates of all 
motors on the lattice towards, or away from the lattice 
tip, respectively. 

The current of particles that leave the system at the 
MT tip depends on depolymerization and detachment 
events and the tip density, 

Jexit(P+) = (5 + /3)p+ ■ (2) 


The tip densities in the IN and the EX phase are obtained 
in a straightforward manner Eng- In the IN phase pb = 
P -, and in the EX phase pb = P+ ■ Bulk and tip currents 
balance in the steady state due to particle conservation, 

JbiPbi P+) = Jexit{p+) > (3) 


and can be solved for the tip density. As results one 
obtains 


in = P(P + ? 7 ~ 1) 

+ p{8 + rf) - 8 - P ’ 


(4) 


and 


P+ X = 1 - 


1 — <5 — 77 


(5) 


In the MC phase the current through the system is de¬ 
termined by the transport capacity of the lattice defined 
from an extremal current principle ISIS], d Pb J b = 0. This 
condition ensues bulk and tip densities in the MC phase 
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FIG. 4. Deviation from the mean-field approximation depending on detachment rates (3, for growth rate p = 0.79, and 
depolymerization rate 5 = 0.2. Panel (a) shows how the data for tip densities deviates from the analytic results for a set of 
system sizes. The deviation from mean-field occurs at the transition between the EX phase and the MC phase, Eqs. fU and 
respectively. Note that the deviation from mean-field result depends on the system size N. (b) Bulk density fluctuations in 
terms of the normalized standard deviation of p\,. An heuristic scaling analysis, where we rescaled the detachment rate as fi'/N, 
shows data collapse for the onset of characteristic density fluctuations, (c) The average current through the system J measured 
from stochastic simulations was evaluated from 10 3 realization of the process, where J = (Q)/At. The number of particles 
that leave the system from the tip through depolymerization or detachment is denoted Q and was measured in time intervals 
At = 10 6 . Scaling is obtained by plotting the data versus . Panel (d) shows the scaled current fluctuations AJ/N of 
the dataset presented in (c), with A J = (( Q 2 ) — {Q) 2 )/At. The data also reveals scaling with fiy/N and supports the scaling 
relation between tip detachment and system size. 


mc _ /3 + # — \/(/3 + 6)(f3 + r)(5 + ??—!)) 

Ph S + r, 

mc = 2/3 + 6(rj + l) + (rj- l)rj - 2y / (/3 + 5)\/3 + ? 7(6 + 77 - 1)] 
P+ {5 + rj) 2 


( 6 ) 

(7) 


To test the validity of the mean-field approach we com¬ 
pare the analytic results for the tip density, p+, with 
stochastic simulation data. For slow growth rates we 
find excellent agreement between the calculations and the 
data, cf. Fig. [3ja) . The transition between EX and IN 
phase is discontinuous as known for TASEP, see Fig. §b). 
For a later comparison with the intermittent regime, we 
also evaluated characteristic density fluctuations across 
the EX/IN transition, see Fig. [3](c). These observations 
are in agreement with the low-density high-density coex¬ 
istence in the classical TASEP [HlH - fTT] , which corresponds 
to the case of S = rj = 0 and p- = /3 in our model. In 


Fig. | density fluctuations show a pronounced peak at 
a critical depolymerization rate 5 C = p_ — /3, which can 
be attributed to the EX/IN transition. Our observation 
of intermittent behavior, suggests that the mean-field ap¬ 
proximation becomes invalid at particular parameter val¬ 
ues. Thereby the detachment rate of motors from the tip, 
/?, plays a critical role. To understand the emergence of 
intermittent behavior, and the eventual break-down of 
the mean-field approximation, we studied the system as 
a function of system size N and /3 in stochastic simula¬ 
tions [J5]. Figure[4](a) shows how the data for the tip den¬ 
sities deviate from the analytic results for a growth 
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rate r) = 0.79 and a depolymerization rate 5 = 0.2. While 
there is good agreement for relatively large motor detach¬ 
ment rates (3 > 0.03, analytic results and the data deviate 
significantly for smaller values of (3. To exclude the pos¬ 
sibility that these deviations from mean-field results are 
due to finite size effects, we also recorded the bulk den¬ 
sity and the particle current through the system as well 
as the fluctuations of these two quantities. As already 
shown above in Fig. [2] the bulk density switches between 
states of high and low density in the intermittent regime. 
Thus, we hypothesize that density and current fluctua¬ 
tions are characteristic for the intermittent regime. To 
test this hypothesis we evaluated the bulk density with 
respect to fluctuations in terms of the normalized stan¬ 
dard deviation y / Var(p b ) = \/( p {;) — (Pb) 2 /(Pb)- The re¬ 
sults are shown in Fig.^b). We found that for all studied 
system sizes, ranging from N = 200 to N = 3200, fluc¬ 
tuations peak at a particular value of (3. This excludes 
the possibility that the observed phenomenon is a finite 
size effect. Furthermore, by plotting the data versus the 
rescaled detachment rate, (3\/N, we found data collapse 
for the onset as well as for the peak of the fluctuations. 
This heuristic analysis ensues the following system size 
dependent law for the onset of intermittent dynamics 

(3 C oc N~ 1 ' 2 . (8) 

The numerical measurement of the average current in the 
system, J , and the current fluctuations A J, confirms the 
above scaling behavior, see Fig.[4jc) and (d). Note, that 
the onset of intermittency also depends on the actual 
values of 5 and 77 , as will be shown later. 


to our work is that in [43] and m the defects are road¬ 
blocks, whereas in our model the kinetic processes at the 
MT tip can function as a defect. 

[4S1 investigated a system with 
defect particles on a constantly growing lattice, where the 
defect property is passed from one particle to the next 
against the direction of transport. This ensues defect 
propagation and interesting DW dynamics in the bulk 
of the lattice. A simple mean-field approach fails to de¬ 
scribe the system, but a careful analysis of the different 
spatial domains that evolve allows an analysis of the de¬ 
fect dynamics. The main differences to our model are, 
the feedback between lattice dynamics and tip occupa¬ 
tion, and that the MT tip site is the only defect in the 
system we study. 

In the following we propose a domain wall theory to 
explore phase diagrams of TASEP systems where tip dy¬ 
namics and bulk dynamics are coupled through shorten¬ 
ing or growth processes. The theory reveals dynamical 
phase transitions in the model and thereby explains why 
simulation data deviate from the mean-field calculations, 
cf. Fig. (4^a). On the level of individual trajectories, 
however, the mean-field results can be verified: The ve¬ 
locities of domain walls can be read out directly from 
kymographs. Our approach generalizes previous work on 
domain wall theory HE EH towards an understanding of 
driven diffusive systems of interacting particles on dy¬ 
namically evolving lattices. 


Pinkoviezky and Gov 


D. Domain wall theory 


C. Formation of shocks 


The appearance and dissolution of motor traffic jams 
at MT tips, provides a mechanism for stochastic switch¬ 
ing between states of growth and shrinking. Recent stud¬ 
ies have investigated TASEP systems with similarly com¬ 
plex shock dynamics, see [43l I46H48] . In the following we 
briefly discuss methods and results of some of those ref¬ 
erences and briefly elaborate on the differences to our 
work. 

|Turci' et al. [23] investigated a system with a defect 
which underlies on/off kinetics on an otherwise static 
lattice with TASEP dynamics. Intermittent density fluc¬ 
tuations are observed, including strong deviations from 
the classic mean-field approach. To improve beyond a 
mean-field approach, the authors employ an intermit¬ 
tent mean-field approach which allows to calculate the 
average current in the system during intermittency. |Sa-| 
hoo et al. 32 investigated a defect which in addition to 
attachment/detachment kinetics also diffuses. In their 
model traffic jams form and dissolve stochastically in the 
bulk of the lattice. The bulk defects considered in the 
above references mm lead to phenomena which are 
similar to our observations at the MT tip. The difference 


Because analytic expressions for the currents, the tip 
densities, and the bulk densities are known for the dif¬ 
ferent phases, we employ a domain wall (DW) theory 
and extremal current principle [38l 149M5T] . In short, DW 
theory determines how shocks and density perturbations 
evolve in the system: The direction in which a shock 
moves and the direction in which a density perturbation 
spreads determine the existence and stability of the dif¬ 
ferent phases [49] . In other words, we are interested in 
the sign of various DW velocities wdWj and the sign the 
collective velocity v co ii which tells about the spreading of 
density perturbations. In the model considered here the 
bulk currents are not independent of the tip density, as in 
the TASEP. Microscopically this means that the occupa¬ 
tion number of the tip n + £ {0,1} influences DW motion 
and thus the phase behavior of the system, see the kymo¬ 
graph in Fig. [5] for example. The figure shows stochastic 
switching in the tip occupation number. Switching be¬ 
tween n+ = 0 and n+ = 1 is indicated by arrows and 
dashed horizontal lines. The solid lines are guides to 
the eye indicating the correlation between tip occupa¬ 
tion and distinct DW velocities in the bulk of the lattice. 
If switching events occur within a relatively short time 
window, multiple DWs coexist in the bulk of the lattice 
as indicated by the numbers (J) and (2) in Fig. [5] The 
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FIG. 5. Kymograph of the model in a regime with strong 
density fluctuations. The data shows how two traffic jams 
nucleate at the MT tip and propagate into the bulk of the 
system before they dissolve. Dashed lines and triangles in¬ 
dicate switching of the tip occupation between n+ = 1 and 
n+ = 0. Shock formation and propagation in the bulk of 
the lattice is highlighted by white solid lines. Parameters are 
P = 0.005, p- = 0.5, 5 = 0.2, 7? = 0.79. 


DW velocity can be analyzed analytically, it is given by 


^yleft _ fright 


^DW = 


pleft _ pright ’ 


(9) 


where left and right denote the densities and currents on 
either side of a DW. The sign of wdw determines whether 
a shock in the system travels to the left (upw < 0 ) or the 
rig ht (upw > 0) . For our purposes the typical procedure 
by Kolomeisky et al. [H] needs to be modified, because 
the currents J on either side of the DW depend explicitly 
on the tip density p+ as illustrated in Fig. [ 6 ] In the 
following we investigate upw between the MC phase and 
the EX phase, since we have seen that this transition is 
involved in the stochastic switching between a growing 
and a shrinking state of the lattice, cf. Fig. [7] For the 
moment, consider the MC phase at the left side of a DW 
and the EX phase at its right. This implies that the tip 
density is in the EX phase p + = p+ x . The DW velocity 
Deft/right then reads 


wmc/ex — 


j(^pr)-wspr) 


p MC 


Pb^ 


( 10 ) 


Note that both currents, J left and J rlght , depend on the 
same tip density, p+ x , but differ with respect to the bulk 
densities p left = p^ G and p rlght = p® x . Vice versa, if the 
MC phase is on the right side of the DW the tip density 
is p+ = p+ c . With the EX phase at the left, the DW 
velocity is 


%X/MC 


wspt 3 )- Apr,pf c ) 


Pb x 


Pb c 


(ii) 


Evaluating the DW velocities in the system by employing 
the current given in Eq. ([!]) and the tip and bulk densities 



FIG. 6. Illustration of possible domain walls (a) and a density 
perturbation (b) in a background density pb- The spike-like 
behavior of the density profile at the lattice tip is also il¬ 
lustrated. Note the particular role of the tip density p+: it 
influences the currents in the system. 


for the individual phases, we are able to construct the 
phase diagram. In the following we provide an attempt 
for a characterization of the arising phases. 


E. Stripe phase 


As suggested from the data shown in the kymograph of 
Fig.0 there is a regime in which the system completely 
switches from the MC phase to the EX phase for small p. 
Because the switching appears as a band of motor parti¬ 
cles across the lattice, we refer to this regime as a stripe 
phase, cf. Fig. [7](a). A reasonable and simple condition 
for the stripe phase is that all three domain walls need 
a negative DW velocity. This is evident from the ky¬ 
mograph of Fig. [7] For ft = 0 this condition cannot be 
fulfilled: 


Dn/ex = —p + 5, 

(12) 

^EX/MC = i 

(13) 

%C/EX = 0 , 

(14) 


where right after the switching event to growth [Eq. (131] 
we assumed p[, lsht = p+ = 0 , and right after the switching 
event to depolymerization [Eq. @] we assumed p^ Ight = 
p + = 1, as suggested from the data. The first two of the 
above equations indeed show a negative DW velocity as 
expected. The MC/EX domain wall, however, does not 
show a negative velocity as observed in the simulations. 
Let us note that for ft = 0 we find a critical growth 
rate r/ s = 1 — S (dashed line in Fig. [ 8 ]) or likewise, a 
critical depolymerization rate = 1 — 77 [S] . To better 
explore the possibility of i'mc/ex < 0 we study the case 
P > 0. We find that for finite values of /3 the critical 
line broadens up into a regime in which Umc/ex < 0 
is possible. In terms of the growth rate, the conditions 
for this regime can be found from i>mc/ex = 0 and that 
%c/ex Fas to be real in the relevant parameter regime. 
It follows that the parameter region of «mc/ex < 0 is 
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FIG. 7. (a) Kymograph of stochastic switching from lattice growth to shrinking and back. Solid lines indicate the velocities 
at which DWs propagate in the bulk of the lattice. For the time points indicated by thin dashed lines the density profiles 
are illustrated in panels (b) and (c) as indicated. Panel (b) shows the DW dynamics before and after the catastrophe event. 
In (c) the DW dynamics before and after the rescue event is illustrated, (d) Illustrates the consecutive phase changes in the 
system which we call “traffic dynamic instability”, indicative of the fact that the system does not settle into its non-equilibrium 
stationary states. Parameters in (a) are /3 = 0.001, = 0.5, 5 = 0.2, rj = 0.79. 
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FIG. 8. Phase diagram for the stripe phase. For f3 = 0 the 
system shows a discontinuous transition between EX and MC 
phase, where stripes emerge right at the critical line (dashed). 
In the case ft > 0 a distinct region in phase space exists in 
which the system robustly forms stripe patterns due to the 
creation of shocks at the MT plus end and a subsequent switch 
to polymerization. 

given by 

V> ^(l-<5+^-4/? + (l-<5) 2 ) , (15) 

l ~ S { l + n-P + s) ' (16) 

The above relations provide the condition for stripe for¬ 
mation in the system based solely on domain wall the¬ 


ory, see shaded are in Fig. [8j Although we find quali¬ 
tative agreement with individual kymographs, compare 
Fig. [TJa) and illustrations in Fig. [7|b), (c) and (d), the 
stochastic nature of the phenomenon makes it hard to 
quantitatively confirm this regime numerically. The ma¬ 
jor difficulty thereby is that the regime of intermittent 
dynamics is not only confined to the stripe phase but 
extends to above and below in phase space. In these 
regimes the system does not reach a non-equilibrium 
steady state either, but is continuously driven between 
different phases in course of time. This view is fortified 
in the following section. 

F. Intermittent phase 

To characterize the non-equilibrium state in which the 
system does not settle into a particular (non-equilibrium) 
stationary state, we distinguish domain walls between 
different phases considering the IN, EX, and MC phase. 
An illustration of the possible domain walls in the system 
is shown in Fig. [bja). In particular, we are interested in 
the direction of DW motion. This quantity determines if 
a transient DW moves to the left or to the right in the 
system. 

Thus, we investigated the sign of iipw for the six pos¬ 
sible combinations of DWs in the system. Our results 
are summarized in Tab. [TJ In the following we discuss 
the cases related to the EX phase, because they are par- 


























ticularly relevant to understand the intermittent regime. 
We begin with the DW between the EX and the MC 
phase. Notably neither of the two phases, MC or EX at 
the right hand side of the system, is stable against the 
other phase on the left hand side of the system. This 
means, that depending on the initial preparation of the 
system, MC/EX or EX/MC, the final state of the system 
is MC or EX, respectively, and thus ergodicity seems to 
be broken. Similar observations were made recently for 
an exclusive queuing model [52] ■ The parameter regime 
where the phenomenon occurs can be determined from 
the conditions Wmc/ex > 0 and Uex/mc > 0 - We find a 
critical growth rate rj c , which reads 


Tj C = 1 S 


1 6 f3 


(17) 


Further, the criterion for the EX/MC phase transition, 
in the absence of particle detachment from the tip, coin¬ 
cides with the change of sign in uex/mc from positive to 
negative: 


77 * = 1-5. (18) 

Consequently for rj > 77 * the MC phase is stable as noted 
previously [ 6 ]. The system behavior for 77 < 77 ° remains 
to be determined. As shown in Tab. U there is also the 
possibility that an IN phase exists in the system below a 
critical growth rate 77 < 77 ^. This finding is in agreement 
with our earlier observation that density and current fluc¬ 
tuations are large in this parameter regime, cf. Fig. [4] 
Here we refrain from a more detailed investigation of the 
IN phase and study the role of fluctuations for the EX 
phase instead, because only if the EX phase is unstable 
intermittency is possible. To this end we analyze the col¬ 
lective velocity in the EX phase. The collective velocity 
is defined as 


recoil — r?pb Jh (pb, P+ ) ■ (19) 

It probes the stability of a small density perturbation in 
a background bulk density, see Fig. [6](b) for an illustra¬ 
tion. Quite generally, the EX phase is characterized by 
negative collective velocity, 7 ;^; < 0, because in a high 
background density a small perturbation moves to the 
left [49] . However, in contrast to typical TASEP sys¬ 
tems, this is not the case in our model. In the regime 
of 77 c < 77 < 77 *, the velocity of a perturbation is positive 
^coii > indicating that density perturbations in the sys¬ 
tem travel to the right. Interestingly the critical growth 
rate which determines tj^ = 0 coincides with the critical 
growth rate as obtained from the DW analysis above, rj c . 

Until now we have relied on the principles of classical 
DW analysis to determine the phase behavior of the sys¬ 
tem. We have learned that between the EX and the MC 
phase, for r] c < 77 < 77 *, DW velocities and the collec¬ 
tive velocity of the EX phase show interesting behavior. 
For 77 > 77 * the system is robustly in the MC phase [ 6 ], 
and for rj < rj c the EX phase is stable in the system. 
In the case if < 77 < 77 *, however, the system does not 


TABLE I. Summary of results for DW motion. rj c re fers to the 
value determined by the EX/MC case given by Eq. |l7| . The 
boxed cases indicate mutually unstable cases of DW motion. 


Domain wall 

v<v c 


?7 C < 77 < 77 * 

. * 

V > V 

MC|EX 

n.a. 




EX|MC 

n.a. 


“T»l 



MC & EX 

MC 

IN|EX 

_AF 


_AF 


EX|IN 



Tt. 

^k 


EX & IN 


EX 

ex & irj“ 

MC|IN 

n.a. 



_iF 

IN|MC 

n.a. 


^k 

^k 




MC 

MC 


a Note this combination is not realistic, because neither the EX 
nor the IN phase exist in this parameter regime [§]. 


settle into a stationary state, but is continuously driven 
between different transient states. So far this behavior 
can be summarized as follows 


77 > 77 * 
rf < 77 < 77 * 
77 < ?7 C 


MC, 

MC & EX, 
EX. 


As discussed in the previous sections, it appears from 
stochastic realizations of the system (Fig. [5]), that per¬ 
turbations at the MT tip render the bulk of the system 
unstable and promote the observed intermittent behav¬ 
ior of the system. This importance of the boundaries 
in driven diffusive systems was recognized by Krug |38| . 
And the model we study here adds dynamic complexity 
to boundary induced phase transitions [38] . Therefore 
we try to include stochastic switching at the MT tip into 
the analysis of the collective velocity as discussed above. 
In the following we show how the collective velocity is af¬ 
fected from switching events at the tip, by assuming that 
the tip density p+ is characterized by only two states, 
given by the tip occupation number, n + = 1 and n+ = 0 . 
Thus we assume that the tip density takes only values 
p+ = 1 and p+ = 0 , which leads to the following equa¬ 
tions for the collective velocity: 


v co n(p+= 1) = 1 - 2p h + 6, ( 20 ) 

Won {p+ = 0 ) = 1 - 2p h - 77 . ( 21 ) 


The above equations illustrate that stochastic switching 
between an empty and an occupied MT tip can promptly 
affect the sign of the collective velocity. This can be seen 
immediately, because depolymerization (rate S) and poly¬ 
merization (rate rf) contribute with different signs to 77 C oii- 
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(a b) 



(depolymerization) S 



0 400 


Site 


FIG. 9. (a) Phase diagram of the model in terms of depolymerization and growth rates for p_ = 0.5 and (5 = 0.05. The shaded 
areas show regimes of intermittent dynamics as calculated in the main text. The topology of the diagram with respect to high 
density (EX), low density (IN) and maximal current density (MC) phases correspond to the case of /3 = 0 [6j. The inset shows 
how bulk density fluctuations are increased in the intermittent regime. A comparison with Fig. [3jc) reveals that fluctuations 
are enhanced by almost a factor of four as compared to the discontinuous EX/IN transition. Data was recorded for system 
size N = 400 and p = 0.4 + <5 along the thick dashed line, (b) Representative kymographs across the intermittent regime. 
Parameter values correspond to those of the filled circles in the inset of panel (a). 


This means that depending on the state of the MT tip 
- occupied or empty - perturbations may either stabilize 
or destabilize the bulk density. This can be illustrated 
in a straight forward manner for the IN phase, when we 
choose pb = p~ =1/2 for example. In this case the sign 
of f co n depends only on the tip density: For p + = 1 and 
Pb = p™ the lattice is in the shrinking state, while par¬ 
ticles on the lattice still travel towards the tip at unit 
velocity. Consequently perturbation travel to the right 
with u co n = 5. This follows directly from Eq. ( |20| and 
is in line with the kymographic data presented above. 
For p+ = 0 the lattice is in the growing state, and par¬ 
ticles on the lattice travel towards the tip at a reduced 
velocity 1 — p. As can be seen from Eq. (21), density 


perturbation travel to the left with negative collective 
velocity u co ii = — 77 , again in agreement with the kymo¬ 
graphic data. Now let us apply this argument for the EX 
phase, where the situation is more intricate. The sign- 
change of the collective velocity, u co n(pj/ x , p+ = 0) < 0 
and u C oii(pjp/p-f- = 1) > 0, is restricted to a particular 
regime rp < p < p c , with 


77 


t=i-a- 


2/3 

1-5 ’ 


( 22 ) 


while the remaining part of the parameter space is not af¬ 
fected. As a consequence the EX phase is destabilized by 
molecular switching events at the MT tip in this regime. 
This means that for growth rates p' < p < p c the molec¬ 
ular noise due to motor occupation at the tip renders the 
EX phase unstable. This mechanism is complementary 


to the above argument for DW velocities. In Fig.[9])a) we 
show the phase diagram with the EX, MC and IN phase 
as well as the different intermittent regimes. In the inset 
we quantified density fluctuations across the intermittent 
region in phase space. The different regimes, p' < p < p c 
(light gray) and p c < p < p* (darker gray), can be identi¬ 
fied reasonably well, given the complexity of the dynam¬ 
ics. As will be shown in the next section it is instrumental 
to employ higher order moments and a direct evaluation 
of bulk density distributions to distinguish between the 
regimes. Figure [9)(b) shows kymographic data of typi¬ 
cal trajectories in the intermittent regime. This direct 
visualization of the process is hard to analyze compu¬ 
tationally, but in the eye of the beholder, domain walls 
and qualitative differences in current fluctuations can be 
readily recognized. 


G. Bistability and collective motion 

To complete the picture of the intermittent regime, 
we performed extensive stochastic simulations to obtain 
the bulk density distribution, see Fig. Uta). We fur¬ 
ther analyze the data using the bimodality parameter 
b = (pl\ + l)//t 4 , where ^4 are the third and fourth 
standardized central moments of the density distribu¬ 
tion. Figure [l0|(b) displays the results. The intermittent 
regime is indicate by light and dark gray regions in the 
plot. At 5\ b increases with increasing 5 and p : it peaks 
at S c and drops to a constant value b = 1/3 for 6 > S*. 
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FIG. 10. (a) Distribution of the bulk density from trajec¬ 
tories as long as 10 s time steps for 0 = 0.05 (thick), (b) 
The bimodality parameter b, supports the role of fluctuations 
during traffic dynamic instability: 6 = 1/3 corresponds to a 
Gaussian distribution and b = 5/9 to a uniform distribution. 
A distribution with 6 > 5/9, as reached for 5 C , can be con¬ 
sidered bimodal. Simulations were conducted at system size 
N = 400 and y = 0.4 + 5. 


The latter value corresponds to a Gaussian distribution. 
The above analysis corroborates our findings for the in¬ 
termittent regime: Stochastic switching between differ¬ 
ent phases ensues large density fluctuations which can 
be recognized in terms of a bistable bulk density distri¬ 
bution. 

Finally, we wish to address the scaling relation we have 
found for density and current fluctuations in the begin¬ 
ning of the paper [Fig. [4] . Although we have discussed 
the origin of fluctuations and intermittent dynamics in 
the system, an explanation for the characteristic scal¬ 
ing 0 C oc TV” 1 / 2 for the onset of density fluctuations has 
remained elusive. In order to gain at least some phe¬ 
nomenological insight, kymographic data was recorded 
for a set of system sizes and 0 = \N ~ 1 / 2 , what approx¬ 
imately corresponds to the onset of the regime of strong 
fluctuation. The results are shown in Fig. El In particu¬ 
lar, we chose the fields of view in a mainly growing state, 
where the system is most likely in the MC phase. A close 
inspection reveals that density fluctuations on the lattice 
do not vanish in a diffusive manner, but rather perform 
erratic zig-zag motion on the lattice, as indicated by the 
symbols. This is opposed to the situation 0 > 0 C , where 
perturbations drift to the left of the system. Similarly 
for 0 < 0 C the shrinking phase becomes more likely and 
perturbations move towards the tip of the lattice. The 
data shown in Fig. |TT] suggests that individual density 
fluctuations are maintained and held in the bulk of the 
system through the combined dynamics of growth and 
shrinking. A detailed analysis however lies beyond the 
scope of this paper. We think that a theory in which 
perturbations were considered as quasi-particles or col¬ 
lective excitations of the system are likely to provide a 
deeper understanding of the TV” 1 / 2 relation. 


IV. DISCUSSION AND CONCLUSION 

In this article we presented a driven lattice gas inspired 
from microtubule depolymerizing motors (kinesin-8), and 
based on the totally asymmetric simple exclusion pro¬ 
cess. The system exhibits boundary-induced intermit¬ 
tent dynamics, with stochastic switching between differ¬ 
ent phases of motor traffic and between lattice growth 
and shrinking. The standard mean-field approach ex¬ 
plains density and current profiles of the system, but 
misses a proper description of the intermittent regime. 
We introduce an extended domain wall theory comple¬ 
mented by an extremal current principle which predicts 
the intermittent regime, in which multiple phases coexist. 
This phase coexistence can not be resolved with time- or 
ensemble-averages, because in such data, it appears the 
system deviates from the usually well-behaved mean-field 
case. However, an interpretation of individual realiza¬ 
tions of the process (like in a single molecule experiment) 
allows to apply the mean-field results to the data. 

Intermittency in driven systems has been reported be¬ 
fore. For example in bidirectional transport |53H55| , and 
in cellular automata for traffic flow [5BJ [57] . The reasons 
for the interesting dynamics in these systems are parti¬ 
cle interactions in bulk of the systems — in contrast to 
the situation considered here, where the system is trig¬ 
gered at the boundary. Quite generally, the presence of 
bottlenecks leads to interesting effects which affect the 
particle densities in driven systems. This includes bot¬ 
tlenecks in bulk Ha m\, at system boundaries [58l - f60] . 
as well as dynamically evolving bottlenecks pRlUT^l . In 
our case we could attribute the intermittent dynamics 
to bulk density fluctuations which arise through molec¬ 
ular noise at the lattice tip, i.e. switching between tip 
occupations n+ = 1 and n+ = 0. The sources of this 
noise is motor detachment and arrival at the tip. In our 
model, the bulk density fluctuations can be attributed 
to a region in phase space in which different phases of 
motor traffic coexist, or in other words, multiple phases 
can be transiently stable for a given set of parameters. 
To identify these dynamical phase transitions and track 
them analytically we neglected attachment and detach¬ 
ment kinetics of the motors, and thus assumed a constant 
density profile. 

In a biological situation where molecular motors inter¬ 
act with MTs, this assumption is valid when the length 
of the filament exceeds a critical length oc w a -1 , where 
w a is the attachment rate of motors to the MT [29]. We 
checked numerically that the intermittent regime also oc¬ 
curs in the case of a constant density profile with explicit 
motor attachment and detachment. This can be under¬ 
stood from a lattice gas point of view, where motor de¬ 
tachment from the tip is equivalent to the creation of a 
“hole” at the tip: A hole can reach the tip by particle 
detachment in bulk and subsequent transport of the hole 
to the tip through depolymerization of the lattice. 

The stochastic switching between growth and shrink¬ 
ing in our model emerges from collective effects of molec- 
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Site 


FIG. 11. At the system size dependent critical detachment rate /3 C (N ), perturbations underlie erratic zig-zag motion in the bulk 
of the system as highlighted by the triangles and vertical dashed lines. Shown are kymographic data for several system sizes 
as indicated. Periods of depolymerization and polymerization alternate in a way that density perturbations are maintained 
within the system. Perturbations seem to perform random walks that are created and annihilated stochastically. Simulation 
parameters were <5 = 0.2, p = 0.79, p_ = 0.5, and, from left to right, /3 C (400) = 0.0125, /3 C ( 800) = 0.008838, /3 C (1600) = 0.00625. 


ular motor traffic. For MT dynamic instability in con¬ 
trast the mechanism of stochastic switching between 
growing and shrinking can be attributed to the nucleotide 
state of the MT lattice [6T]. A comparison by numbers 
could be thought of within the mathematical framework 
provided by Dogteronr and Leibler [62] , where the param¬ 
eters are growth and shrinking speeds z>+, V-, and the fre¬ 
quencies of catastrophe /_) and rescue /_| Our model 

relates microscopically to these macroscopic parameters. 
The speed of the MT tip is a function of the motor density 
at the MT tip and given by v(p + ) = p(l — p + ) — Sp + SB]- 
Depending on whether the system is in a high density 
phase or a low density phase v(p + ) can be negative or 
positive. For the switching frequencies the mapping is 
only possible for the rescue frequency. It can be di¬ 
rectly attributed to the parameter of particle detachment 
from the tip /_ oc /3. Catastrophe in contrast is ini¬ 

tiated by particle arrival at the tip and initiation of a 
traffic jam. Thereby current fluctuations are important 
and more elaborate techniques are necessary to calculate 
these [Hj. 

A future challenge is to understand the interplay be¬ 
tween the MT lattice, MT tips, and MT associated 
proteins. The characteristics of each of these parts 


contribute significantly to MT dynamics, but the re¬ 
lations between them are still obscure; necessary in¬ 
formation is available: protein localization mechanisms 
at MT tips [63], enzymatic functions of tip related en¬ 
zymes [641165] , and dynamic information about the MT 
tip structure mm- 

The present model constitutes an example how en¬ 
zymes may influence MT dynamics. The mechanism we 
found is complementary to known mechanisms of MT 
regulation, but may contribute to MT regulation in a 
similar way as does MT dynamic instability [5]. 
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